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Abstract 

Use of the Smagorinsky eddy-viscosity formulation and related schemes for 
sub-grid-scale parameterization of large eddy simulation models requires speci- 
fication of a single length scale, earlier related by Lilly to the scale of filtering 
and/or numerical resolution. An anisotropic integration of the Kolmogoroff en- 
strophy spectrum allows generalization of that relationship to anisotropic res- 
olution. It is found that the Deardorff assumption is reasonably accurate for 
small anisotropies and can be simply improved for larger values. 

1. Introduction 

Numerical integration of the time-dependent equations of fluid dynamics to 
simulate evolution of the largest scales of motion has been applied for about 30 
years, initially in the field of weather prediction. Early workers were concerned 
over what to do about the scales of motion too small to be resolved by the 
computer, but the problem was often confused with the numerical errors intro- 
duced by finite difference algorithms. The technique that Smagorinsky (1963) 
introduced grew out of an empirical variable eddy viscosity method due, I be- 
lieve to R. Richtmyer (no reference available), for smoothing simulated shock 
wave calculations. Smagorinsky and his associates recognized, however, that 
Richtmyer’s variable viscosity was also consistent with the notion of a universal 
equilibrium range of turbulence. I pursued this point and used the Kolmogoroff 
inertial sub-range hypothesis to quantitatively relate the length scale required 
in the Smagorinsky formulation to the resolved scale (Lilly, 1966, 1967). I as- 
sumed isotropic resolution, although numerical simulations then and now often 
are carried out with anisotropic resolution. The question of how to generalize 
my expression for this purpose has engaged some attention. Deardorff (1970) 
assumed that the length scale is proportional to the cube root of the product 
of the resolution scales in the three directions. Alternatives have also been pro- 
posed, most of them involving the effects of a nearby boundary. Piomelli, et al 
(1987) review the subject and conclude that no real consensus has yet emerged, 
though Deardorff ’s formulation is widely used. 
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The purpose of this note is to aid in resolving this issue. I extend my earlier 
calculation, oriented now toward spectral and pseudo-spectral integration tech- 
niques, to allow resolution to vary with direction. I still assume, however, that 
the small scale limit of the simulation is in the Kolmogoroff inertial sub-range in 
all directions. This is a serious limitation, as the most common reason for apply- 
ing anisotropic resolution is an expectation of anisotropic and/or inhomogeneous 
turbulence, typically in regions close to a boundary. Thus the scale anisotropy 
problem is usually compounded with a boundary layer problem and with the 
likelihood that a classical inertial sub-range may not exist, at least in the scales 
resolved by the simulation. Piomelli, et al (1987) discuss and test several addi- 
tional parameterizations intended to improve the boundary layer resolution, but 
all are modifications of an isotropic formulation related to Deardorff’s. The re- 
sults to be presented here are offered to test and perhaps improve the Deardorff 
assumption. 


2. The problem and the solution 

Consider the problem of integrating the Navier-Stokes equation under condi- 
tions of incompressibility and constant density, so that 


dui /dt + (9(uiti/) / dxj -f dp/dxi = ud 2 Ui jdx 2 (l) 


and 


du{/dxi= 0 , ( 2 ) 

where p is pressure divided by density, and other terminology is conventional. A 
low-pass spatial filter is applied to these equations, in recognition of the limited 
resolution available to any real discretization technique. In the early develop- 
ment era this was typically a uniform average over a box surrounding a spatial 
grid point. Leonard (1974) introduced the use of Gaussian filters, a technique 
which is widely applied in engineering fluid dynamics simulations. In many 
current research simulations the effective filter is the spectral cut-off of a finite 
Fourier transform. In any case it is assumed to be a linear operator and is here 
designated by angle brackets <>. Because it is linear it may be exchanged with 
the differential operators in (l) and (2) and applied to the flow variables directly. 
This allows the linear terms of the equations to be integrated as if they were 
unfiltered. The filtered momentum flux product, < Ujtiy >, cannot, however, be 
resolved directly from the filtered velocities, and must be subject to a creative 
parameterization . 

I assume the decomposition of the momentum flux product applied by Sma- 
gorinsky, Lilly, and Deardorff, i.e. 
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< U,Uy >—< Ui >< U y > —Tij + 

Tij = — < U,Uy > + < Ui >< Uj > —SijE 

E = (<u 2 > - < Ui > 2 )/2 (3) 

with Tij designated as the subgrid scale (SGS) stress, E the SGS kinetic energy, 
and 6ij is the Kronecker delta. The filtered equations of motion and continuity 
are now written as 

d < Ui > /dt + 3(< Ui >< Uj >)/dxj + dn/dxi 

= drij/dxj + ud 2 < Ui > /dx 2 (4) 

d < Ui > /dxi — 0, (5) 

where i r =< p > +E. 

The Smagorinsky parameterization for the SGS stress is of the eddy viscosity 
type, given by 

nj = K(Sij + Sji), ( 6 ) 

where S^y — d < Ui > /dxj, and 

K = A 2 5, (7) 

where S 2 = S»y(S,y + Sji). The length scale A was assumed by Smagorinsky to 
be proportional to the grid resolution. Some investigators have replaced S by 
lo, the square root of the enstrophy, i.e. to 2 = (3u, / <3xy — 3uy/3x t ) 2 / 2. Not 
much difference is found in practice, corresponding to the fact that the volume 
integrals of S 2 and w 2 are identical and they have the same Fourier spectral 
amplitudes. 

I evaluated the length scale A by assuming the validity of inertial cascade the- 
ory and neglecting intermittency effects. If the resolved kinetic energy equation 
is derived by multiplying (4) by < u,- >, the contribution from the eddy stress 
term is 

< Ui > drij/dxj = d(< Ui > Tij)/dxj — KS 2 

= d(< Ui > Tij)/dxj - A 2 S 3 (8) 

The first term on the rhs is variable in sign and vanishes or is small in the 
volume average, while the second term is always negative. It is assumed that 
the loss of energy from the resolved scales, A 2 S 3 , is the same as the viscous 
dissipation of energy, e. These are equated after evaluating 5 2 from spectral 
integration, i. e. 

S 2 = [ m k 2 E(k)dk 


( 9 ) 
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where E(k) is the spectral amplitude of kinetic energy, defined as an integral 
along the surfaces of spheres in spectral space. The wavenumber k m must cor- 
respond to the resolution limit, which is the cut-off wavenumber for spectral 
discretization. It is assumed that at and near k m the energy spectrum is given 
by the inertial range form, i.e. 

E(k) = ae 2 / 3 fc- 6 / 3 , (10) 

where a is the Kolmogoroff constant, about 1.5, and e is the viscous dissipation 
of kinetic energy. Substitution of (10) into (9) leads to the result of integration 

S 2 = (3/4 )ae 2 / 3 fc m 4 ' 3 (ll) 

The rate of energy removal from the resolved flow is evaluated by substituting 

(ll) into the last term of (8) so that 

A 2 S 3 = (3a/4) 3 / 2 A 2 ek m 2 (12) 

Equating this expression to dissipation allows evaluation of the length scale as 

A = S~ 3 /V/2 = ( 4 /3a) 3 / < ‘fe m - 1 (13) 

One may now set k m = constant/ A, where A is the grid interval, assumed 
to be isotropic. The minimum value for the constant, assuming the spectral 
wavelengths are bounded by the Nyquist limit, is two. When pseudo-spectral 
integration techniques are applied, in order to avoid aliasing the constant must 
be no less than three, which increases A by 50%. 

I now wish to drop the assumption of resolution isotropy. For simplicity 
I assume, however, axisymmetry, i.e. that the limits of resolution in two of 
the three dimensions are the same. The typical situation involves increased 
resolution close to a boundary, either through a variable grid length or use 
of Chebyshev polynomials. Within the axisymmetric framework one also may 
consider the effects of decreased resolution in one direction. In the conclusion 
section I argue that the results for different resolution in all three directions can 
be determined adequately from the axisymmetric case. 

Results are obtained by writing the enstrophy spectrum function in three- 
dimensional spectral space and integrating it over prolate or oblate spheroids, 
that is figures obtained by rotating an ellipse about its long or short axis. This 
may not correspond accurately to the resolution limits in a model using Fourier 
modes in two dimensions and Chebyshev polynomials or finite differences in the 
third, but it is more readily compared with the isotropic case. Similar results 
can be obtained by assuming a cylindrical volume in wavenumber space, but 
they differ only slightly and not in any qualitatively important manner. 
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If the three-dimensional spectrum function is isotropic, as assumed, it is 
E(k) / 4nk 2 , and the 3-d enstrophy or S 2 spectrum is E(k)/4ir. The wavenumber 
spheroid is assumed to be bounded by the equation for an ellipse, i.e. 

{kh/khm) 2 + ( k x /k zm ) 2 = 1 

where kh is the wavenumber in the axisymmetric direction, notationally assumed 
to be horizontal, and k x is that in the vertical direction, with kh m and k xm the 
horizontal and vertical spectral limits. In place of (9), the evaluation of S 2 is 
accomplished by integration over the spheroid as follows: 

S 2 = / / E{k)k h dk h dk x (14) 

Jo Jo 

The limit on the first integral is khi = ^ m (l — k x 2 /k xm 2 ) 1 / 2 - Note that if 
E{k) were replaced by 4n in (14), the integral would give the volume of the 
wavenumber spheroid, 47rfcj, m 2 fc zm . With the energy spectrum given by (10), 
with k 2 = kh 2 + k 2 , integration of (14) yields 

S 2 = (3/4 )ae 2/5 k hm 

where r = k hm /k xm and y = r 10 / 9 

The particular form for S is chosen so that if it is substituted into the first 
equality of (13), with y = 1 and the wavenumbers assumed to be inversely pro- 
portional to the grid spacings Ax, Ay, A z, one obtains the expression assumed 
by Deardorff, i.e. A ~ (AxAyAz) 1 / 3 . More generally, upon substituting S into 
(13) one obtains 

X = A/>y -3 / 4 , Xd = (4/3a) 3 ^ 4 fcfc m 2 ^ 3 k xm ~ 1 ^, (16) 

with Xjo equivalent to Deardorff ’s expression. 

The integral for y in (15) can be reduced to the forms 

y = r */ 9 {l-r 2 )- 1 ' 2 r^cosx^dx.tan*™ = (1 - r 2 ) 1 / a /r,for r < 1 (17a) 
Jo 

y = r 4 / 9 (r 2 - l)" 1 / 2 I*™ (cos x)~ 2 t z dx, sin x m = (r 2 - l) x / 2 /r,for r > 1 (176) 

Jo 

For the extreme cases r = 0, oo, the integration limits x m become 7r/2 in both 
expressions, and the integrals are evaluated in terms of tabulated Gamma func- 
tions. Thus for small r, S 2 ~ khm*/ 5 and A ~ khm 1 , while for large r, 
S 2 ~ kh m 1 / 3 *. 

zm and A ~ khm 1 '*k xm 3 ^ 4 - A somewhat more revealing in- 
terpretation can be obtained by differentiating A with respect to r. This yields 
the following: 


8/9 fc^ 4/S y(r) (15) 

[ [r 2 + (1 - r 2 )x 2 ] -6/6 dx. 
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— 1/3 for r C 1 (18a) 

r\~ 1 d\/dr = —(3 r/4y)dy/dr = 0 for r = 1 (186) 

5/12 for r > 1 (18c) 

d 2 \/d(lnr) 2 = (4/27)A r> at r = 1. (18 d) 

These show that the Deardorff expression is valid in the vicinity of r = 1, but 
that with increasing anisotropy in either direction A becomes larger than A d- 
Fig. 1 is a plot of A/ A £>. The data for the curve labelled “Length scale” were 
calculated to two-three digit accuracy by using the (I8d) for .05 < r < 20 
and matching that with an expansion of the expressions in (l7a,b) around their 
limiting values. The straight lines are plots of the asymptotic forms, that is 
y~ 3 / 4 for the y's given by (17a) for r C 1 and (17b) for r> 1. 

An accurate approximation for .02 < r < 50 is 

A/A 0 = r( 2 / 27 ) 1/, +r-( 2 / 27 ) 1/J -l (19) 

This is obtained as a solution to (l8d), but with the rhs replaced by (2/27) (A -f- 
A.d). This relation is also plotted in Fig. 1, and is imperceptibly different from 
the “Length scale” except at the most extreme anisotropies. 

3. Conclusions 

The results of the above calculation indicate that the Deardorff formulation 
is accurate to within 20% for .2 < r < 5. For anisotropies greater than that the 
length scale should be increased, with (19) sufficiently accurate for anisotropies 
less than a factor of 50 either way. Such large anisotropies are unlikely to produce 
accurate simulations in any case, and probably need to be enhanced by more 
sophisticated parameterizations than the Smagorinsky formulation. 

Upon consideration of the nature of these solutions, it seems evident that 
they can be extended to the case where the resolution is different in all three 
directions. Since equation (18) or (19) is symmetric in Inr, it doesn’t seem to 
matter much whether two wavenumber limits are larger and one smaller, or one 
larger and two smaller. Apparently the results are only sensitive to the largest 
ratio of the length scales, so presumably the existence of a third dimension with 
an intermediate resolution would have little effect. 
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FIGURE 1. The turbulence length scale, normalized by DeardorfTs (1970) 
assumption, and determined by solution of (13) and (15), for various ratios of 
the limiting wavenumbers, k^m/kzm- The straight lines area symptotic forms 
from (l7a,b) and the curve marked with + signs is an approximation given by 
(19). 



